Tasks Complete
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(stringr)
library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
gfr <- read.csv("eGFR.csv")
weight <- read.csv("weight.csv")
weight$Visit <- as.numeric(str_extract(weight$Study.visit, "[0-9]"))
weight$Id <- as.numeric(sapply(str_split(weight$Id, "-"), function(v) v[3]))
weight$Pregnant[which(weight$Pregnant == "")] <- NA
weight$Date <- mdy(weight$Date)
gfr <- left_join(gfr, weight)
## Joining, by = c("Id", "Visit")
gfr$Monday <- wday(gfr$Date, label = TRUE) == "Mon"
gfr <- gfr %>%
mutate(Season = case_when(Season == "Rainy season" ~ "Rainy", Season == "Dry season" ~
"Dry")) %>%
mutate(SeasonRainy = if_else(Season == "Rainy", 1, 0))
gfr <- gfr %>%
mutate(Pregnant = if_else(is.na(Pregnant), "No", Pregnant)) %>%
mutate(PregnantYes = if_else(Pregnant == "Yes", 1, 0)) %>%
select(Id, Visit, Sex, Community, Age, eGFR_cre, Months, Height..meters., Weight..kg.,
BMI, Season, SeasonRainy, Pregnant, PregnantYes, Date, Monday)
head(gfr)
## Id Visit Sex Community Age eGFR_cre Months Height..meters. Weight..kg. BMI
## 1 1 2 Male Tololar 1 22 123.8020 6 1.78 80.3 25.34
## 2 1 6 Male Tololar 1 25 118.3518 36 1.78 78.9 24.90
## 3 1 4 Male Tololar 1 23 116.7217 18 1.78 82.7 26.10
## 4 1 5 Male Tololar 1 24 122.0748 24 1.78 77.6 24.49
## 5 1 7 Male Tololar 1 27 122.9960 48 1.78 79.4 25.06
## 6 1 3 Male Tololar 1 23 128.7687 12 1.78 72.6 22.91
## Season SeasonRainy Pregnant PregnantYes Date Monday
## 1 Dry 0 No 0 2015-05-27 FALSE
## 2 Rainy 1 No 0 2017-10-23 TRUE
## 3 Dry 0 No 0 2016-06-03 FALSE
## 4 Rainy 1 No 0 2016-10-17 TRUE
## 5 Rainy 1 No 0 2018-10-22 TRUE
## 6 Rainy 1 No 0 2015-10-29 FALSE
Analytic Subset (n = 301)
New variables
gfr_sub <- gfr %>%
group_by(Id) %>%
mutate(Count = n()) %>%
ungroup() %>%
filter(Count > 3) %>%
select(-Count)
gfr_sub <- gfr_sub %>%
filter(Visit == 1 & Pregnant == "No") %>%
group_by(Id) %>%
mutate(Baseline = mean(eGFR_cre)) %>%
ungroup() %>%
distinct(Baseline, Id) %>%
right_join(gfr_sub) %>%
filter(!is.na(Baseline))
## Joining, by = "Id"
gfr_sub %>%
distinct(Id, Baseline) %>%
ggplot(aes(x = Baseline)) + geom_histogram() + geom_vline(xintercept = 90) +
theme_classic()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
gfr_sub2 <- gfr_sub %>%
# filter(Baseline > 60) %>%
mutate(egfrChg = eGFR_cre - Baseline) %>%
mutate(egfrPer = (eGFR_cre - Baseline)/Baseline) %>%
select(Id, Months, egfrChg, egfrPer, eGFR_cre, Baseline, Sex, Pregnant, PregnantYes,
Season, SeasonRainy, Age, Monday, Weight..kg., Visit) %>%
arrange(Id, Months)
Observation: Individuals have slightly lower eGFR for their normal on Monday than other days.
# Monday vs. Other Days
gfr_sub2 %>%
group_by(Id) %>%
mutate(zegfr = scale(eGFR_cre)) %>%
ggplot(aes(x = Monday, y = zegfr)) + geom_boxplot() + ylab("eGFR Subject Z Score") +
theme_classic()
Observation: Individuals have higher eGFR in comparison to their normal when pregnant.
gfr_sub2 %>%
filter(Sex == "Female") %>%
group_by(Id) %>%
mutate(zegfr = scale(eGFR_cre)) %>%
ggplot(aes(x = Pregnant, y = zegfr)) + geom_boxplot() + ylab("eGFR Subject Z Score") +
theme_classic()
Observation: No relationship between weight and eGFR.
gfr_sub2 %>%
group_by(Id) %>%
mutate(zegfr = scale(eGFR_cre)) %>%
ggplot(aes(x = Weight..kg., y = zegfr, color = Sex)) + geom_point() + xlab("Weight (kg)") +
ylab("eGFR Subject Z Score") + theme_classic()
Observation: No relationship between weight and eGFR.
gfr_sub2 %>%
group_by(Id) %>%
mutate(zweight = scale(Weight..kg.)) %>%
mutate(zegfr = scale(eGFR_cre)) %>%
ggplot(aes(x = zweight, y = zegfr, color = Sex)) + geom_point() + xlab("Weight Subject Z Score") +
ylab("eGFR Subject Z Score") + theme_classic()
Observation: Individuals have higher eGFR in comparison to their normal during Dry Season
gfr_sub2 %>%
group_by(Id) %>%
mutate(zegfr = scale(eGFR_cre)) %>%
ggplot(aes(x = Season, y = zegfr)) + geom_boxplot() + ylab("eGFR Subject Z Score") +
theme_classic()
gfr_sub2 %>%
ggplot(aes(x = Months, y = egfrPer, group = Id)) + geom_line() + theme_classic()
gfr_sub2 %>%
ggplot(aes(x = Months, y = egfrChg, group = Id)) + geom_line() + theme_classic()
gfr_sub2 %>%
filter(Visit == 1) %>%
nrow()
## [1] 303
gfr_sub2 %>%
filter(Visit == 1) %>%
count(Sex)
## # A tibble: 2 x 2
## Sex n
## <chr> <int>
## 1 Female 75
## 2 Male 228
gfr_sub2 %>%
filter(Visit == 1) %>%
pull(Age) %>%
summary()
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.00 21.00 24.00 23.71 27.00 30.00
gfr_sub2 <- gfr_sub2 %>%
mutate(Years = Months/12) %>%
mutate(eGFRCat = factor(case_when(eGFR_cre > 90 ~ "Normal", eGFR_cre <= 90 ~
"Low"))) %>%
mutate(eGFRPerCat = factor(case_when(egfrPer > -0.1 ~ "Normal", egfrPer <= -0.1 ~
"Low")))
gfr_sub2 %>%
filter(Id == 172) %>%
ggplot(aes(x = Years, y = eGFR_cre)) + geom_line() + theme_classic()
gfr_sub2 %>%
filter(Id == 121) %>%
ggplot(aes(x = Years, y = eGFR_cre)) + geom_line() + theme_classic()
gfr_sub2 %>%
filter(Id == 212) %>%
ggplot(aes(x = Years, y = eGFR_cre)) + geom_line() + theme_classic()
gfr_sub2 <- gfr_sub2 %>%
filter(!(Id %in% c(172, 121)))
#Model Estimation
Outcome Assumptions
State Assumptions
Transition Assumptions
Transition intensity from state r to s is the e instantaneous risk of moving from state r to state s, modeled as
\[q_{rs}(z) = q_{rs}^{(0)}\exp(\beta^T_{rs}z)\]
and \(q_{rr} = -\sum_{s\not= r}q_{rs}\).
Constraints
library("msm")
hmodel1 <- list(hmmNorm(mean = 0, sd = 5), hmmNorm(mean = -30, sd = 8))
mod1 <- msm(egfrChg ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0,
0.01), c(0.01, 0)), hmodel = hmodel1, covariates = ~eGFRCat, initprobs = c(0.95,
0.05), control = list(reltol = 1e-25, maxit = 10000), fixedpars = c(7, 8))
Fitted Model
Baseline Transition Intensities = \(q_{rs}^{(0)}\)
Hazard ratios = \(\exp(\beta_{rs})\)
Parameters
State 1: Healthy
State 2: Unhealthy
mod1
##
## Call:
## msm(formula = egfrChg ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0, 0.01), c(0.01, 0)), hmodel = hmodel1, covariates = ~eGFRCat, initprobs = c(0.95, 0.05), fixedpars = c(7, 8), control = list(reltol = 1e-25, maxit = 10000))
##
## Maximum likelihood estimates
## Baselines are with covariates set to their means
##
## Transition intensities with hazard ratios for each covariate
## Baseline eGFRCatNormal
## State 1 - State 1 -0.09873 (-0.12646,-0.07708)
## State 1 - State 2 0.09873 ( 0.07708, 0.12646) 0.8167 (0.4068,1.640)
## State 2 - State 1 0.20529 ( 0.11602, 0.36324) 2.4352 (0.8296,7.148)
## State 2 - State 2 -0.20529 (-0.36324,-0.11602)
##
## Hidden Markov model, 2 states
## State 1 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean -3.033914 -3.405789 -2.662039
## sd 7.603481 7.325144 7.892394
##
## State 2 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean -30 NA NA
## sd 8 NA NA
##
##
## -2 * log-likelihood: 16940.57
## [Note, to obtain old print format, use "printold.msm"]
Estimated probability of going from state r to s in the period of t years for given covariate values
library(purrr)
library(tidyr)
library(tibble)
dat <- crossing(t = 0:5, A = 0:1)
tranprobs <- function(x, t, A) {
as.vector(pmatrix.msm(x = x, t = t, covariates = list(eGFRCatNormal = A)))
}
dat %>%
select(t, A) %>%
pmap(.f = tranprobs, x = mod1) %>%
enframe() %>%
mutate(transition = list(c("State 1 to State 1", "State 2 to State 1", "State 1 to State 2",
"State 2 to State 2"))) %>%
bind_cols(dat %>%
select(t, A)) %>%
unnest() %>%
rename(prob = value) %>%
mutate(eGFR = case_when(A == 1 ~ "Normal eGFR", TRUE ~ "Low eGFR")) %>%
# mutate(Sex = case_when(B == 1 ~ 'Male', TRUE ~'Female')) %>%
ggplot(aes(x = t, y = prob, color = transition, group = transition)) + geom_line() +
facet_grid(. ~ eGFR) + ylim(c(0, 1)) + labs(y = "Probabilitiy of Transition in Time Period",
x = "Time Period (Years)", color = "Transition") + theme_classic()
## Warning: `cols` is now required when using unnest().
## Please use `cols = c(value, transition)`
hazard.msm(mod1) #hazard ratios with CI
## $eGFRCatNormal
## HR L U
## State 1 - State 2 0.8167364 0.4068397 1.639609
## State 2 - State 1 2.4351988 0.8296341 7.147962
Outcome Assumptions
State Assumptions
Transition Assumptions
Transition intensity from state r to s is the e instantaneous risk of moving from state r to state s, modeled as
\[q_{rs}(z) = q_{rs}^{(0)}\exp(\beta^T_{rs}z)\]
and \(q_{rr} = -\sum_{s\not= r}q_{rs}\).
Constraints
library("msm")
hmodel2 <- list(hmmNorm(mean = 0, sd = 0.1), hmmNorm(mean = -0.3, sd = 0.1))
mod2 <- msm(egfrPer ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0,
0.01), c(0.01, 0)), hmodel = hmodel2, covariates = ~eGFRCat, initprobs = c(0.95,
0.05), control = list(reltol = 1e-25, maxit = 10000), fixedpars = c(7, 8))
Fitted Model
Baseline Transition Intensities = \(q_{rs}^{(0)}\)
Hazard ratios = \(\exp(\beta_{rs})\)
Parameters
State 1: Healthy
State 2: Unhealthy
mod2
##
## Call:
## msm(formula = egfrPer ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0, 0.01), c(0.01, 0)), hmodel = hmodel2, covariates = ~eGFRCat, initprobs = c(0.95, 0.05), fixedpars = c(7, 8), control = list(reltol = 1e-25, maxit = 10000))
##
## Maximum likelihood estimates
## Baselines are with covariates set to their means
##
## Transition intensities with hazard ratios for each covariate
## Baseline eGFRCatNormal
## State 1 - State 1 -0.07623 (-0.09973,-0.05827)
## State 1 - State 2 0.07623 ( 0.05827, 0.09973) 0.2596 (0.1434,0.4699)
## State 2 - State 1 0.11207 ( 0.03983, 0.31539) 1.1354 (0.2437,5.2887)
## State 2 - State 2 -0.11207 (-0.31539,-0.03983)
##
## Hidden Markov model, 2 states
## State 1 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean -0.02542415 -0.02902379 -0.02182450
## sd 0.07318805 0.07051590 0.07596146
##
## State 2 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean -0.3 NA NA
## sd 0.1 NA NA
##
##
## -2 * log-likelihood: -4150.321
## [Note, to obtain old print format, use "printold.msm"]
Estimate probability of going from state r to s in the period of t years for given covariate values
library(purrr)
library(tidyr)
library(tibble)
dat <- crossing(t = 0:5, A = 0:1)
tranprobs <- function(x, t, A) {
as.vector(pmatrix.msm(x = x, t = t, covariates = list(eGFRCatNormal = A)))
}
dat %>%
select(t, A) %>%
pmap(.f = tranprobs, x = mod2) %>%
enframe() %>%
mutate(transition = list(c("State 1 to State 1", "State 2 to State 1", "State 1 to State 2",
"State 2 to State 2"))) %>%
bind_cols(dat %>%
select(t, A)) %>%
unnest() %>%
rename(prob = value) %>%
mutate(eGFR = case_when(A == 1 ~ "Normal eGFR", TRUE ~ "Low eGFR")) %>%
# mutate(Sex = case_when(B == 1 ~ 'Male', TRUE ~'Female')) %>%
ggplot(aes(x = t, y = prob, color = transition, group = transition)) + geom_line() +
facet_grid(. ~ eGFR) + ylim(c(0, 1)) + labs(y = "Probabilitiy of Transition in Time Period",
x = "Time Period (Years)", color = "Transition") + theme_classic()
## Warning: `cols` is now required when using unnest().
## Please use `cols = c(value, transition)`
hazard.msm(mod2) #hazard ratios with CI
## $eGFRCatNormal
## HR L U
## State 1 - State 2 0.2595931 0.1434145 0.4698868
## State 2 - State 1 1.1353537 0.2437311 5.2887307
Outcome Assumptions
State Assumptions
Transition Assumptions
Transition intensity from state r to s is the e instantaneous risk of moving from state r to state s, modeled as
\[q_{rs}(z) = q_{rs}^{(0)}\exp(\beta^T_{rs}z)\]
and \(q_{rr} = -\sum_{s\not= r}q_{rs}\).
Constraints
library("msm")
hmodel3 <- list(hmmNorm(mean = 120, sd = 12), hmmNorm(mean = 60, 15))
mod3 <- msm(eGFR_cre ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0,
0.01), c(0.01, 0)), hmodel = hmodel3, hcovariates = list(~PregnantYes + SeasonRainy,
~PregnantYes + SeasonRainy), hconstraint = list(PregnantYes = c(1, 1)), initprobs = c(0.95,
0.05), control = list(reltol = 1e-25, maxit = 10000), center = FALSE) #,fixedpars=c(5,6)
Fitted Model
Baseline Transition Intensities = \(q_{rs}^{(0)}\)
Hazard ratios = \(\exp(\beta_{rs})\)
Parameters
State 1: Healthy
State 2: Unhealthy
mod3
##
## Call:
## msm(formula = eGFR_cre ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0, 0.01), c(0.01, 0)), hmodel = hmodel3, hcovariates = list(~PregnantYes + SeasonRainy, ~PregnantYes + SeasonRainy), hconstraint = list(PregnantYes = c(1, 1)), initprobs = c(0.95, 0.05), center = FALSE, control = list(reltol = 1e-25, maxit = 10000))
##
## Maximum likelihood estimates
## Baselines are with covariates set to 0
##
## Transition intensities
## Baseline
## State 1 - State 1 -0.07321 (-0.099392,-0.053929)
## State 1 - State 2 0.07321 ( 0.053929, 0.099392)
## State 2 - State 1 0.00815 ( 0.001338, 0.049651)
## State 2 - State 2 -0.00815 (-0.049651,-0.001338)
##
## Hidden Markov model, 2 states
## State 1 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean 121.911144 121.004627 122.8176621
## sd 9.066679 8.608588 9.5491469
## PregnantYes 13.418100 7.933667 18.9025323
## SeasonRainy -1.615515 -2.679358 -0.5516715
##
## State 2 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean 78.489625 73.9837877 82.995461
## sd 23.514248 22.1634408 24.947384
## PregnantYes 13.418100 7.9336669 18.902532
## SeasonRainy 4.506299 -0.3606806 9.373279
##
##
## -2 * log-likelihood: 17816.77
## [Note, to obtain old print format, use "printold.msm"]
Estimate probability of going from state r to s in the period of t months for given covariate values
library(purrr)
library(tidyr)
library(tibble)
dat <- crossing(t = 0:3, A = 0:1)
tranprobs <- function(x, t, A) {
as.vector(pmatrix.msm(x = x, t = t, covariates = list(eGFRPerCatNormal = A)))
}
dat %>%
select(t, A) %>%
pmap(.f = tranprobs, x = mod3) %>%
enframe() %>%
mutate(transition = list(c("State 1 to State 1", "State 2 to State 1", "State 1 to State 2",
"State 2 to State 2"))) %>%
bind_cols(dat %>%
select(t, A)) %>%
unnest() %>%
rename(prob = value) %>%
mutate(eGFR = case_when(A == 1 ~ "Normal eGFR Change", TRUE ~ "High eGFR Decline")) %>%
# mutate(Sex = case_when(B == 1 ~ 'Male', TRUE ~'Female')) %>%
ggplot(aes(x = t, y = prob, color = transition, group = transition)) + geom_line() +
facet_grid(. ~ eGFR) + ylim(c(0, 1)) + labs(y = "Probabilitiy of Transition in Time Period",
x = "Time Period (Years)", color = "Transition") + theme_classic()
hazard.msm(mod3) #hazard ratios with CI
## [1] "No covariates on transition intensities"
Outcome Assumptions
State Assumptions
Transition Assumptions
Transition intensity from state r to s is the e instantaneous risk of moving from state r to state s, modeled as
\[q_{rs}(z) = q_{rs}^{(0)}\exp(\beta^T_{rs}z)\]
and \(q_{rr} = -\sum_{s\not= r}q_{rs}\).
Constraints
library("msm")
hmodel4 <- list(hmmNorm(mean = 120, sd = 12), hmmNorm(mean = 60, sd = 15))
mod4 <- msm(eGFR_cre ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0,
0.01), c(0.01, 0)), hmodel = hmodel4, hcovariates = list(~PregnantYes + SeasonRainy,
~PregnantYes + SeasonRainy), hconstraint = list(PregnantYes = c(1, 1)), covariates = ~eGFRPerCat,
initprobs = c(0.95, 0.05), control = list(reltol = 1e-25, maxit = 10000), center = FALSE)
Fitted Model
Baseline Transition Intensities = \(q_{rs}^{(0)}\)
Hazard ratios = \(\exp(\beta_{rs})\)
Parameters
State 1: Healthy
State 2: Unhealthy
mod4
##
## Call:
## msm(formula = eGFR_cre ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0, 0.01), c(0.01, 0)), hmodel = hmodel4, covariates = ~eGFRPerCat, hcovariates = list(~PregnantYes + SeasonRainy, ~PregnantYes + SeasonRainy), hconstraint = list(PregnantYes = c(1, 1)), initprobs = c(0.95, 0.05), center = FALSE, control = list(reltol = 1e-25, maxit = 10000))
##
## Maximum likelihood estimates
## Baselines are with covariates set to 0
##
## Transition intensities with hazard ratios for each covariate
## Baseline eGFRPerCatNormal
## State 1 - State 1 -0.15471 (-0.297238,-0.080521)
## State 1 - State 2 0.15471 ( 0.080521, 0.297238) 3.956e-01 (0.179,0.8742)
## State 2 - State 1 0.02105 ( 0.004787, 0.092545) 5.689e-07 (0.000, Inf)
## State 2 - State 2 -0.02105 (-0.092545,-0.004787)
##
## Hidden Markov model, 2 states
## State 1 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean 121.814771 120.899092 122.7304490
## sd 9.139974 8.688090 9.6153608
## PregnantYes 13.512675 7.985020 19.0403289
## SeasonRainy -1.609828 -2.669647 -0.5500089
##
## State 2 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean 77.860355 73.3378509 82.382860
## sd 23.444099 22.0848117 24.887048
## PregnantYes 13.512675 7.9850202 19.040329
## SeasonRainy 4.711851 -0.1885915 9.612294
##
##
## -2 * log-likelihood: 17810.79
## [Note, to obtain old print format, use "printold.msm"]
Estimate probability of going from state r to s in the period of t months for given covariate values
library(purrr)
library(tidyr)
library(tibble)
dat <- crossing(t = seq(0, 3, by = 0.5), A = 0:1)
tranprobs <- function(x, t, A) {
as.vector(pmatrix.msm(x = x, t = t, covariates = list(eGFRPerCatNormal = A)))
}
dat %>%
select(t, A) %>%
pmap(.f = tranprobs, x = mod4) %>%
enframe() %>%
mutate(transition = list(c("State 1 to State 1", "State 2 to State 1", "State 1 to State 2",
"State 2 to State 2"))) %>%
bind_cols(dat %>%
select(t, A)) %>%
unnest() %>%
rename(prob = value) %>%
mutate(eGFR = case_when(A == 1 ~ "Normal eGFR Change", TRUE ~ "High eGFR Decline")) %>%
# mutate(Sex = case_when(B == 1 ~ 'Male', TRUE ~'Female')) %>%
ggplot(aes(x = t, y = prob, color = transition, group = transition)) + geom_line() +
facet_grid(. ~ eGFR) + ylim(c(0, 1)) + labs(y = "Probabilitiy of Transition in Time Period",
x = "Time Period (Years)", color = "Transition") + theme_classic()
hazard.msm(mod4) #hazard ratios with CI
## $eGFRPerCatNormal
## HR L U
## State 1 - State 2 3.955948e-01 0.1790104 0.8742245
## State 2 - State 1 5.688522e-07 0.0000000 Inf
Outcome Assumptions
State Assumptions
Transition intensity from state r to s is the e instantaneous risk of moving from state r to state s, modeled as
\[q_{rs}(z) = q_{rs}^{(0)}\exp(\beta^T_{rs}z)\]
and \(q_{rr} = -\sum_{s\not= r}q_{rs}\).
Constraints
library("msm")
hmodel5 <- list(hmmNorm(mean = 120, sd = 12), hmmNorm(mean = 70, sd = 15))
mod5 <- msm(eGFR_cre ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0,
0.01), c(0.01, 0)), hmodel = hmodel5, hcovariates = list(~PregnantYes + SeasonRainy,
~PregnantYes + SeasonRainy), hconstraint = list(PregnantYes = c(1, 1)), covariates = ~eGFRPerCat,
initprobs = c(0.95, 0.05), control = list(reltol = 1e-25, maxit = 10000), fixedpars = c(7,
8), center = FALSE)
Fitted Model
Baseline Transition Intensities = \(q_{rs}^{(0)}\)
Hazard ratios = \(\exp(\beta_{rs})\)
Parameters
State 1: Healthy
State 2: Unhealthy
mod5
##
## Call:
## msm(formula = eGFR_cre ~ Years, subject = Id, data = gfr_sub2, qmatrix = rbind(c(0, 0.01), c(0.01, 0)), hmodel = hmodel5, covariates = ~eGFRPerCat, hcovariates = list(~PregnantYes + SeasonRainy, ~PregnantYes + SeasonRainy), hconstraint = list(PregnantYes = c(1, 1)), initprobs = c(0.95, 0.05), fixedpars = c(7, 8), center = FALSE, control = list(reltol = 1e-25, maxit = 10000))
##
## Maximum likelihood estimates
## Baselines are with covariates set to 0
##
## Transition intensities with hazard ratios for each covariate
## Baseline eGFRPerCatNormal
## State 1 - State 1 -0.13315 (-0.22648,-0.07828)
## State 1 - State 2 0.13315 ( 0.07828, 0.22648) 3.268e-01 (0.1673,0.6385)
## State 2 - State 1 0.05451 ( 0.02411, 0.12322) 4.627e-06 (0.0000, Inf)
## State 2 - State 2 -0.05451 (-0.12322,-0.02411)
##
## Hidden Markov model, 2 states
## State 1 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean 120.983101 120.025993 121.9402087
## sd 9.923965 9.517420 10.3478762
## PregnantYes 14.368859 8.473357 20.2643603
## SeasonRainy -1.856296 -2.952367 -0.7602244
##
## State 2 - normal distribution
## Parameters:
## Estimate LCL UCL
## mean 70.000000 NA NA
## sd 15.000000 NA NA
## PregnantYes 14.368859 8.473357 20.264360
## SeasonRainy 6.809362 4.649229 8.969496
##
##
## -2 * log-likelihood: 18069.84
## [Note, to obtain old print format, use "printold.msm"]
Estimate probability of going from state r to s in the period of t months for given covariate values
library(purrr)
library(tidyr)
library(tibble)
dat <- crossing(t = seq(0, 3, by = 0.5), A = 0:1)
tranprobs <- function(x, t, A) {
as.vector(pmatrix.msm(x = x, t = t, covariates = list(eGFRPerCatNormal = A)))
}
dat %>%
select(t, A) %>%
pmap(.f = tranprobs, x = mod5) %>%
enframe() %>%
mutate(transition = list(c("State 1 to State 1", "State 2 to State 1", "State 1 to State 2",
"State 2 to State 2"))) %>%
bind_cols(dat %>%
select(t, A)) %>%
unnest() %>%
rename(prob = value) %>%
mutate(eGFR = case_when(A == 1 ~ "Normal eGFR Change", TRUE ~ "High eGFR Decline")) %>%
# mutate(Sex = case_when(B == 1 ~ 'Male', TRUE ~'Female')) %>%
ggplot(aes(x = t, y = prob, color = transition, group = transition)) + geom_line() +
facet_grid(. ~ eGFR) + ylim(c(0, 1)) + labs(y = "Probabilitiy of Transition in Time Period",
x = "Time Period (Years)", color = "Transition") + theme_classic()
hazard.msm(mod5) #hazard ratios with CI
## $eGFRPerCatNormal
## HR L U
## State 1 - State 2 3.268474e-01 0.167305 0.6385298
## State 2 - State 1 4.627078e-06 0.000000 Inf
Comparison of predicted States for each time across all individuals at https://bheggeseth.shinyapps.io/KidneyModelExploration/
tmp1 <- viterbi.msm(mod1)
tmp1$State <- apply(tmp1$pstate, 1, which.max)
tmp2 <- viterbi.msm(mod2)
tmp2$State <- apply(tmp2$pstate, 1, which.max)
tmp3 <- viterbi.msm(mod3)
tmp3$State <- apply(tmp3$pstate, 1, which.max)
tmp4 <- viterbi.msm(mod4)
tmp4$State <- apply(tmp4$pstate, 1, which.max)
tmp5 <- viterbi.msm(mod5)
tmp5$State <- apply(tmp5$pstate, 1, which.max)
# tmp1 <- viterbi.own(mod1,gfr_sub2) tmp2 <- viterbi.own(mod2,gfr_sub2) tmp3 <-
# viterbi.own(mod3,gfr_sub2) tmp4 <- viterbi.own(mod4,gfr_sub2) tmp5 <-
# viterbi.own(mod5,gfr_sub2)
## `summarise()` has grouped output by 'Id', 'Baseline'. You can override using the `.groups` argument.
## Joining, by = "Id"
Visualizations of subject-specific trajectorires
vizTraj <- function(SBJ, tmp, title) {
tmp$Certainty <- apply(tmp$pstate, 1, max)
tmp %>%
mutate(State = factor(fitted)) %>%
filter(subject %in% SBJ) %>%
left_join(gfr_sub2 %>%
select(Id, Years, eGFR_cre, Pregnant, Sex), by = c(subject = "Id", time = "Years")) %>%
ggplot(aes(x = time, y = eGFR_cre, group = subject, color = State)) + geom_line() +
geom_hline(yintercept = c(60, 90), color = "lightgrey", linetype = 3) + geom_point(aes(shape = Pregnant,
alpha = Certainty)) + facet_wrap(Sex ~ subject) + scale_alpha(limits = c(0,
1)) + xlab("Years of Follow Up") + ylab("eGFR") + ylim(c(0, 150)) + scale_color_discrete("State",
limits = factor(2:1), labels = c("Unhealthy", "Healthy")) + ggtitle(title) +
theme_classic() + theme(legend.position = "top") + guides(alpha = FALSE)
}
i = 1
vizTraj(1:30, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 1-30")) %>%
print()
vizTraj(31:60, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 31-60")) %>%
print()
vizTraj(61:90, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 61-90")) %>%
print()
vizTraj(91:120, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 91-120")) %>%
print()
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(121:150, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 121-150")) %>%
print()
vizTraj(151:180, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 151-180")) %>%
print()
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(181:210, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 181-210")) %>%
print()
vizTraj(211:230, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 211-230")) %>%
print()
Visualizations of subject-specific trajectorires
vizTraj <- function(SBJ, tmp, title) {
tmp$Certainty <- apply(tmp$pstate, 1, max)
tmp %>%
mutate(State = factor(fitted)) %>%
filter(subject %in% SBJ) %>%
left_join(gfr_sub2 %>%
select(Id, Years, eGFR_cre, Pregnant, Sex), by = c(subject = "Id", time = "Years")) %>%
ggplot(aes(x = time, y = eGFR_cre, group = subject, color = State)) + geom_line() +
geom_hline(yintercept = c(60, 90), color = "lightgrey", linetype = 3) + geom_point(aes(shape = Pregnant,
alpha = Certainty)) + facet_wrap(Sex ~ subject) + scale_alpha(limits = c(0,
1)) + xlab("Years of Follow Up") + ylab("eGFR") + ylim(c(0, 150)) + scale_color_discrete("State",
limits = factor(2:1), labels = c("Unhealthy", "Healthy")) + ggtitle(title) +
theme_classic() + theme(legend.position = "top") + guides(alpha = FALSE)
}
i = 2
vizTraj(1:30, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 1-30")) %>%
print()
vizTraj(31:60, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 31-60")) %>%
print()
vizTraj(61:90, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 61-90")) %>%
print()
vizTraj(91:120, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 91-120")) %>%
print()
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(121:150, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 121-150")) %>%
print()
vizTraj(151:180, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 151-180")) %>%
print()
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(181:210, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 181-210")) %>%
print()
vizTraj(211:230, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 211-230")) %>%
print()
Visualizations of subject-specific trajectorires
vizTraj <- function(SBJ, tmp, title) {
tmp$Certainty <- apply(tmp$pstate, 1, max)
tmp %>%
mutate(State = factor(fitted)) %>%
filter(subject %in% SBJ) %>%
left_join(gfr_sub2 %>%
select(Id, Years, eGFR_cre, Pregnant, Sex), by = c(subject = "Id", time = "Years")) %>%
ggplot(aes(x = time, y = eGFR_cre, group = subject, color = State)) + geom_line() +
geom_hline(yintercept = c(60, 90), color = "lightgrey", linetype = 3) + geom_point(aes(shape = Pregnant,
alpha = Certainty)) + facet_wrap(Sex ~ subject) + scale_alpha(limits = c(0,
1)) + xlab("Years of Follow Up") + ylab("eGFR") + ylim(c(0, 150)) + scale_color_discrete("State",
limits = factor(2:1), labels = c("Unhealthy", "Healthy")) + ggtitle(title) +
theme_classic() + theme(legend.position = "top") + guides(alpha = FALSE)
}
i = 3
vizTraj(1:30, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 1-30")) %>%
print()
vizTraj(31:60, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 31-60")) %>%
print()
vizTraj(61:90, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 61-90")) %>%
print()
vizTraj(91:120, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 91-120")) %>%
print()
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(121:150, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 121-150")) %>%
print()
vizTraj(151:180, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 151-180")) %>%
print()
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(181:210, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 181-210")) %>%
print()
vizTraj(211:230, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 211-230")) %>%
print()
Visualizations of subject-specific trajectorires
vizTraj <- function(SBJ, tmp, title) {
tmp$Certainty <- apply(tmp$pstate, 1, max)
tmp %>%
mutate(State = factor(fitted)) %>%
filter(subject %in% SBJ) %>%
left_join(gfr_sub2 %>%
select(Id, Years, eGFR_cre, Pregnant, Sex), by = c(subject = "Id", time = "Years")) %>%
ggplot(aes(x = time, y = eGFR_cre, group = subject, color = State)) + geom_line() +
geom_hline(yintercept = c(60, 90), color = "lightgrey", linetype = 3) + geom_point(aes(shape = Pregnant,
alpha = Certainty)) + facet_wrap(Sex ~ subject) + scale_alpha(limits = c(0,
1)) + xlab("Years of Follow Up") + ylab("eGFR") + ylim(c(0, 150)) + scale_color_discrete("State",
limits = factor(2:1), labels = c("Unhealthy", "Healthy")) + ggtitle(title) +
theme_classic() + theme(legend.position = "top") + guides(alpha = FALSE)
}
i = 4
vizTraj(1:30, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 1-30")) %>%
print()
vizTraj(31:60, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 31-60")) %>%
print()
vizTraj(61:90, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 61-90")) %>%
print()
vizTraj(91:120, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 91-120")) %>%
print()
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(121:150, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 121-150")) %>%
print()
vizTraj(151:180, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 151-180")) %>%
print()
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(181:210, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 181-210")) %>%
print()
vizTraj(211:230, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 211-230")) %>%
print()
Visualizations of subject-specific trajectorires
vizTraj <- function(SBJ, tmp, title) {
tmp$Certainty <- apply(tmp$pstate, 1, max)
tmp %>%
mutate(State = factor(fitted)) %>%
filter(subject %in% SBJ) %>%
left_join(gfr_sub2 %>%
select(Id, Years, eGFR_cre, Pregnant, Sex), by = c(subject = "Id", time = "Years")) %>%
ggplot(aes(x = time, y = eGFR_cre, group = subject, color = State)) + geom_line() +
geom_hline(yintercept = c(60, 90), color = "lightgrey", linetype = 3) + geom_point(aes(shape = Pregnant,
alpha = Certainty)) + facet_wrap(Sex ~ subject) + scale_alpha(limits = c(0,
1)) + xlab("Years of Follow Up") + ylab("eGFR") + ylim(c(0, 150)) + scale_color_discrete("State",
limits = factor(2:1), labels = c("Unhealthy", "Healthy")) + ggtitle(title) +
theme_classic() + theme(legend.position = "top") + guides(alpha = FALSE)
}
i = 5
vizTraj(1:30, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 1-30")) %>%
print()
vizTraj(31:60, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 31-60")) %>%
print()
vizTraj(61:90, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 61-90")) %>%
print()
vizTraj(91:120, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 91-120")) %>%
print()
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(121:150, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 121-150")) %>%
print()
vizTraj(151:180, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 151-180")) %>%
print()
## Warning: Removed 1 rows containing missing values (geom_point).
vizTraj(181:210, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 181-210")) %>%
print()
vizTraj(211:230, get(paste0("tmp", i)), paste0("Model ", i, ": Subjects 211-230")) %>%
print()